slope <- map_dfr(rigal_trends,
  ~.x %>%
  select(siteid, linear_slope),
.id = "response"
)
slope_df <- slope %>%
  pivot_wider(names_from = "response", values_from = "linear_slope")

0.1 Neutral difference

1 Rare species hypothesis

1.1 hypothesis

  • Knowing that high evenness and turnover trends happens where there are no trends in species richness

  • Turnover linked to shorted lived species

  • To highly abundant but stochastic species

  • Community that contains species more abundant globally have more turnover.

  • Community that contains species less occurent globally have more turnover. -> CWM(rarity)

2 Characterise variable species

  • Most abundant when present = average abundance
  • Rarely present = proportion occurence in sampling
sp_mean_abun <- filtered_dataset$measurement %>%
  left_join(filtered_dataset$site_quali, by = "siteid") %>%
  left_join(filtered_dataset$location, by = "siteid") %>%
  group_by(unitabundance, ecoregion, species) %>%
  summarise(summ =
    list(
      enframe(
        c(
          summary_distribution(abundance, na.rm = TRUE)
        )
      )
    )
    ) %>%
  unnest(cols = summ) %>%
  filter(name == c("mean")) %>%
  pivot_wider(names_from = "name", values_from = "value")
#> `summarise()` has grouped output by 'unitabundance', 'ecoregion'. You can
#> override using the `.groups` argument.
species_fq <- filtered_dataset$measurement %>%
  left_join(filtered_dataset$site_quali, by = "siteid") %>%
  left_join(filtered_dataset$location, by = "siteid") %>%
  group_by(unitabundance, ecoregion, species) %>%
  summarise(fq = n()) %>%
  ungroup()
#> `summarise()` has grouped output by 'unitabundance', 'ecoregion'. You can
#> override using the `.groups` argument.

op_fq <- filtered_dataset$measurement %>%
  left_join(filtered_dataset$site_quali, by = "siteid") %>%
  left_join(filtered_dataset$location, by = "siteid") %>%
  group_by(unitabundance, ecoregion) %>%
  summarise(n_op = n()) %>%
  ungroup()
#> `summarise()` has grouped output by 'unitabundance'. You can override using the
#> `.groups` argument.

species_fq_op <- species_fq %>%
  left_join(op_fq, by = c("unitabundance", "ecoregion")) %>%
  mutate(occurence = fq / n_op) %>%
  select(-fq, -n_op) %>%
  ungroup()
ti <- ecdf(filtered_dataset$measurement$abundance)

get_proba_values_from_dist <- function(
  x = NULL,
  emp_pdf = NULL) {
  emp_pdf(x)
}
emp_pdf <- sp_mean_abun %>%
  group_by(unitabundance, ecoregion) %>%
  summarise(emp_pdf = list(ecdf(mean))) %>%
  ungroup()
#> `summarise()` has grouped output by 'unitabundance'. You can override using the
#> `.groups` argument.
rarity_species <- sp_mean_abun %>%
  ungroup() %>%
  left_join(emp_pdf, by = c("unitabundance", "ecoregion")) %>%
  mutate(
    prob_abun = map2_dbl(mean, emp_pdf,
      ~get_proba_values_from_dist(x = .x, emp_pdf = .y))
    ) %>%
  left_join(species_fq_op, by = c("unitabundance", "ecoregion", "species")) %>%
  select(-emp_pdf)
rarity_species %>%
  ggplot(aes(y = prob_abun, x = log(occurence))) +
  geom_point()


rarity_species %>%
  ggplot(aes(x = prob_abun)) +
  geom_histogram()
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.


rarity_species %>%
  ggplot(aes(x = log(occurence))) +
  geom_histogram()
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

2.1 Occurrence

library(ecocom)
cwm_rarity <- filtered_dataset$measurement %>%
  left_join(select(filtered_dataset$site_quali, siteid, unitabundance), by = "siteid") %>%
  left_join(select(filtered_dataset$location, siteid, ecoregion), by = "siteid") %>%
  left_join(
    rarity_species,
    by = c("ecoregion", "unitabundance", "species")
    ) %>%
  group_by(op_id) %>%
  summarise(
    cwm_occurence = calc_cw_mean(
      trait = log(occurence),
      weight = abundance),
    cwm_prob_abun = calc_cw_mean(
      trait = prob_abun,
      weight = abundance
  )
    )
site_median_cwm_rarity <- cwm_rarity %>%
  left_join(select(filtered_dataset$measurement, siteid, op_id), by = "op_id") %>%
group_by(siteid) %>%
  summarise(across(where(is.double), median))
slope_df_rarity <- slope_df %>%
  left_join(site_median_cwm_rarity, by = "siteid")
slope_df_rarity %>%
  ggplot(aes(
      x = cwm_occurence,
      y = jaccard,
      color = cwm_prob_abun
      )) +
  scale_color_viridis() +
  geom_point() +
  geom_smooth(method = "gam")
#> `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'


slope_df_rarity %>%
  ggplot(aes(
      x = cwm_prob_abun,
      y = jaccard,
      color = cwm_occurence
      )) +
  scale_color_viridis() +
  geom_point() +
  geom_smooth(method = "gam")
#> `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'

tar_load(com_mat_site)

stab_sync <- com_mat_site %>%
  mutate(
    richness = purrr::map(mat, compute_sp_nb_from_com_mat),
    richness_med = purrr::map_dbl(richness, median),
    avg_sp = purrr::map(mat, colMeans),
    cov_mat = purrr::map(mat, cov),
    var_sp = purrr::map(cov_mat, diag),
    synchrony = purrr::map_dbl(cov_mat, compute_synchrony),
    cv_sp = purrr::map2_dbl(avg_sp, var_sp, compute_avg_cv_sp),
    cv_com = compute_cv_com(synchrony = synchrony, cv_sp = cv_sp),
    cv_classic = purrr::map2_dbl(cov_mat, mat, function(variance, abundance) {
      sqrt(sum(variance)) / mean(rowSums(abundance))
      })
  ) %>%
select(-starts_with("mat"))
stab_sync %>%
  ggplot(aes(y = 1/cv_com, x = richness_med)) +
  geom_point() +
  geom_smooth(method = "lm")
#> `geom_smooth()` using formula 'y ~ x'
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).


stab_sync %>%
  ggplot(aes(y = synchrony, x = log(richness_med))) +
  geom_point() +
  geom_smooth(method = "lm")
#> `geom_smooth()` using formula 'y ~ x'


stab_sync %>%
  ggplot(aes(y = cv_sp, x = richness_med)) +
  geom_point() +
  geom_smooth(method = "lm")
#> `geom_smooth()` using formula 'y ~ x'

slope_df_stab <- slope_df %>%
  left_join(stab_sync, by = "siteid")
slope_df_stab %>%
  ggplot(aes(y = jaccard, x = cv_sp)) +
  geom_point()

slope_df_stab %>%
  ggplot(aes(y = jaccard, x = synchrony)) +
  geom_point()

slope_df_stab %>%
  ggplot(aes(y = jaccard, x = 1 / cv_com)) +
  geom_point()

3 PCAAAAAAAAAAAAAAAAAAAAAA

library(ade4)
library(factoextra)
res.pca <- dudi.pca(select(slope_df, -siteid),
  scannf = FALSE,   # Hide scree plot
  nf = 5            # Number of components kept in the results
)
fviz_eig(res.pca)

p_pca <- map(list(c(1,2), c(2, 3), c(1, 3)),
  ~fviz_pca_var(res.pca,
    axes = .x,
    col.var = "contrib", # Color by contributions to the PC
    gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
    repel = TRUE     # Avoid text overlapping
    ) +
  theme(legend.position = "none") + labs(title = "")
)

plot_grid(plotlist = p_pca, ncol = 1)
#> Warning: ggrepel: 8 unlabeled data points (too many overlaps). Consider
#> increasing max.overlaps

median_site <- analysis_dataset %>%
  select(
    siteid,
    starts_with("chao_"),
    evenness,
    total_abundance
  ) %>%
  group_by(siteid) %>%
  summarise(across(where(is.double), median, na.rm =TRUE)) %>%
  rename_with(~paste0("med_", .x), where(is.double)) %>%
  left_join(site_median_cwm_rarity, by = "siteid") %>%
  left_join(select(stab_sync, siteid, cv_com, synchrony, cv_sp), by = "siteid") %>%
  left_join(select(slope_df, siteid, hillebrand, jaccard, total), by = "siteid")
pca_turnover <- dudi.pca(
  select(na.omit(median_site), -siteid),
  scannf = FALSE,   # Hide scree plot
  nf = 5            # Number of components kept in the results
)
fviz_eig(pca_turnover)

p_1_2 <- fviz_pca_var(pca_turnover,
  col.var = "contrib", # Color by contributions to the PC
  gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
  repel = TRUE     # Avoid text overlapping
)
p_1_2

3.1 Jaccard distrribution vs classification_station

slope_classification_station <- map_dfr(
  rigal_trends,
  ~.x %>%
    select(siteid, linear_slope, shape_class, direction),
  .id = "response"
)
get_shape_class_x_trends_y <- function(dataset = NULL, x = NULL, y = NULL,
  type = "shape_class") {

  trends <- dataset %>%
    filter(response == y) %>%
    select(all_of(c("siteid", "linear_slope")))

  shape <- dataset %>%
    filter(response == x) %>%
    select(all_of(c("siteid", type)))

  output <- trends %>%
    left_join(shape, by = "siteid")
  return(output)
}
boxplot_trends_shape <- function(
dataset = NULL, x = NULL, y = NULL, type = "shape_class",
geom_fun = geom_boxplot
) {
  give.n <- function(x){
    return(c(y = median(x)*1.25, label = length(x))) 
    # experiment with the multiplier to find the perfect position
}

  get_shape_class_x_trends_y(
    dataset = slope_classification_station,
    x = x,
    y = y,
    type = type 
    ) %>%
  ggplot(aes_string(x = type, y = "linear_slope")) +
  geom_fun() +
  labs(y = paste0("Trends of ", y), x = paste0("Classification of ", x)) +
  stat_summary(fun.data = give.n, geom = "text", fun = median, 
    position = position_dodge(width = 0.75)) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

}
get_shape_class_x_trends_y(
  dataset = slope_classification_station,
  x = "total_abundance",
  y = "jaccard"
) %>%
  ggplot(aes(x = shape_class, y = linear_slope)) +
  geom_boxplot()

boxplot_trends_shape(
  dataset = slope_classification_station,
  x = "total_abundance",
  y = "jaccard",
  type = "direction"
)

boxplot_trends_shape(
  dataset = slope_classification_station,
  x = "species_nb",
  y = "total",
  type = "direction"
)

3.1.1 Hillebrand

3.1.1.1 Species richness

p <- map(c("hillebrand", "total"),
  ~boxplot_trends_shape( dataset =
    slope_classification_station,
  x = "species_nb",
  y = .x,
  type = "shape_class"
))

plot_grid(plotlist = p)

p <- map(c("hillebrand", "total"),
  ~boxplot_trends_shape( dataset =
    slope_classification_station,
  x = "species_nb",
  y = .x,
  type = "direction"
))

plot_grid(plotlist = p)

3.1.1.2 Total abundance

p <- map(c("hillebrand", "total"),
  ~boxplot_trends_shape( dataset =
    slope_classification_station,
  x = "total_abundance",
  y = .x,
  type = "shape_class"
))

plot_grid(plotlist = p)

p <- map(c("hillebrand", "total"),
  ~boxplot_trends_shape( dataset =
    slope_classification_station,
  x = "total_abundance",
  y = .x,
  type = "direction"
))

plot_grid(plotlist = p)

3.1.2 Baselga

3.1.2.1 Species richness

p <- map(c("nestedness", "turnover", "jaccard_dis"),
  ~boxplot_trends_shape( dataset =
    slope_classification_station,
  x = "species_nb",
  y = .x,
  type = "shape_class"
))

plot_grid(plotlist = p)

p <- map(c("nestedness", "turnover", "jaccard_dis"),
  ~boxplot_trends_shape( dataset =
    slope_classification_station,
  x = "species_nb",
  y = .x,
  type = "direction"
))

plot_grid(plotlist = p)

3.1.2.2 Total abundance

p <- map(c("nestedness", "turnover", "jaccard_dis"),
  ~boxplot_trends_shape( dataset =
    slope_classification_station,
  x = "total_abundance",
  y = .x,
  type = "direction"
))

plot_grid(plotlist = p)

p <- map(c("nestedness", "turnover", "jaccard_dis"),
  ~boxplot_trends_shape( dataset =
    slope_classification_station,
  x = "total_abundance",
  y = .x,
  type = "shape_class"
))

plot_grid(plotlist = p)

3.1.3 Jaccard

3.1.3.1 Species richness

boxplot_trends_shape(
  dataset = slope_classification_station,
  x = "species_nb",
  y = "jaccard",
  type = "direction"
)

boxplot_trends_shape(
  dataset = slope_classification_station,
  x = "species_nb",
  y = "jaccard",
  type = "shape_class"
)

3.1.3.2 Total abundance

boxplot_trends_shape(
  dataset = slope_classification_station,
  x = "total_abundance",
  y = "jaccard",
  type = "direction"
)

boxplot_trends_shape(
  dataset = slope_classification_station,
  x = "total_abundance",
  y = "jaccard",
  type = "shape_class"
)

p <- map(c("nestedness", "turnover", "jaccard_dis"),
  ~boxplot_trends_shape( dataset =
    slope_classification_station,
  x = "species_nb",
  y = .x,
  type = "shape_class"
))

plot_grid(plotlist = p)

3.1.4 Total turnover, appearance, disappearance

var_y <- c("total", "appearance", "disappearance")
p <- map(var_y,
  ~boxplot_trends_shape( dataset =
    slope_classification_station,
  x = "species_nb",
  y = .x,
  type = "direction"
))

plot_grid(plotlist = p)

p <- map(var_y,
  ~boxplot_trends_shape( dataset =
    slope_classification_station,
  x = "species_nb",
  y = .x,
  type = "shape_class"
))

plot_grid(plotlist = p)

3.1.4.1 Total abundance

p <- map(var_y,
  ~boxplot_trends_shape( dataset =
    slope_classification_station,
  x = "total_abundance",
  y = .x,
  type = "direction"
))

plot_grid(plotlist = p)

p <- map(var_y,
  ~boxplot_trends_shape( dataset =
    slope_classification_station,
  x = "total_abundance",
  y = .x,
  type = "shape_class"
))

plot_grid(plotlist = p)

4 In stable richness

stable_rich_siteid <- slope_classification_station %>%
  filter(response == "species_nb", direction == "stable") %>%
  .[["siteid"]]
high_jaccard_slope <-
  slope_df %>%
  filter(siteid %in% stable_rich_siteid) %>%
  select(siteid, total_abundance, species_nb, jaccard) %>%
  arrange(desc(jaccard))

site_desc_jaccard_slope <-
  high_jaccard_slope %>%
  .[["siteid"]]
p_high_jaccard_slope <- map(site_desc_jaccard_slope[1:20],
  ~plot_community_data(
    analysis_dataset %>%
    filter(siteid == .x) %>%
    select(siteid, year, jaccard),
    y = "jaccard", x = "year",
  smoothing_method = "lm"))
plot_grid(plotlist = p_high_jaccard_slope, ncol = 3)

p_pop_high_jaccard <-
  map(site_desc_jaccard_slope[1:20],
    ~plot_temporal_population(
      filtered_dataset$measurement %>%
        filter(siteid == .x)
      ) +
    theme(legend.position = "none")
    )
#> Note: Using an external vector in selections is ambiguous.
#> i Use `all_of(species_var)` instead of `species_var` to silence this message.
#> i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
#> This message is displayed once per session.
#> Note: Using an external vector in selections is ambiguous.
#> i Use `all_of(y_var)` instead of `y_var` to silence this message.
#> i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
#> This message is displayed once per session.
#> Note: Using an external vector in selections is ambiguous.
#> i Use `all_of(species)` instead of `species` to silence this message.
#> i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
#> This message is displayed once per session.
plot_grid(plotlist = p_pop_high_jaccard, ncol = 3)
#> Warning: Removed 16 rows containing missing values (position_stack).
#> Warning: Removed 10 rows containing missing values (position_stack).
#> Warning: Removed 15 rows containing missing values (position_stack).
#> Warning: Removed 224 rows containing missing values (position_stack).
#> Warning: Removed 342 rows containing missing values (position_stack).
#> Warning: Removed 24 rows containing missing values (position_stack).
#> Warning: Removed 169 rows containing missing values (position_stack).
#> Warning: Removed 8 rows containing missing values (position_stack).
#> Warning: Removed 50 rows containing missing values (position_stack).
#> Warning: Removed 30 rows containing missing values (position_stack).
#> Warning: Removed 36 rows containing missing values (position_stack).
#> Warning: Removed 24 rows containing missing values (position_stack).
#> Warning: Removed 18 rows containing missing values (position_stack).
#> Warning: Removed 36 rows containing missing values (position_stack).
#> Warning: Removed 2 rows containing missing values (position_stack).
#> Warning: Removed 45 rows containing missing values (position_stack).

p_turnover_high_jaccard_slope <- map(site_desc_jaccard_slope[1:20],
  ~plot_community_data(
    analysis_dataset %>%
    filter(siteid == .x) %>%
    select(siteid, year, turnover),
    y = "turnover", x = "year",
  smoothing_method = "lm"))
plot_grid(plotlist = p_turnover_high_jaccard_slope, ncol = 3)

p_nestedness_high_jaccard_slope <- map(site_desc_jaccard_slope[1:20],
  ~plot_community_data(
    analysis_dataset %>%
    filter(siteid == .x) %>%
    select(siteid, year, nestedness),
    y = "nestedness", x = "year",
  smoothing_method = "lm"))
plot_grid(plotlist = p_nestedness_high_jaccard_slope, ncol = 3)

p_appearance_high_jaccard_slope <- map(site_desc_jaccard_slope[1:20],
  ~plot_community_data(
    analysis_dataset %>%
    filter(siteid == .x) %>%
    select(siteid, year, appearance),
    y = "appearance", x = "year",
  smoothing_method = "lm"))
plot_grid(plotlist = p_appearance_high_jaccard_slope, ncol = 3)

p_disappearance_high_jaccard_slope <- map(site_desc_jaccard_slope[1:20],
  ~plot_community_data(
    analysis_dataset %>%
    filter(siteid == .x) %>%
    select(siteid, year, disappearance),
    y = "disappearance", x = "year",
  smoothing_method = "lm"))
plot_grid(plotlist = p_disappearance_high_jaccard_slope, ncol = 3)

diff_baseline <- analysis_dataset %>%
    group_by(siteid) %>%
    arrange(year) %>%
    summarise(
      diff_baseline = year[2] - year[1],
      span = year[length(year)] - year[1] + 1
    )
slope_df %>%
  left_join(diff_baseline, by = "siteid") %>%
  ggplot(aes(x = diff_baseline, y = jaccard)) +
  geom_point() +
  geom_smooth(method = "gam")
#> `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'

slope_df %>%
  left_join(diff_baseline, by = "siteid") %>%
  ggplot(aes(x = span, y = jaccard)) +
  geom_point() +
  geom_smooth(method = "gam")
#> `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'

  • Lets see communities that becomes more dissimilar:
p_jaccard_low_jaccard <-
  map(site_desc_jaccard_slope[(length(site_desc_jaccard_slope)-20):length(site_desc_jaccard_slope)],
    ~plot_community_data(
      analysis_dataset %>%
        filter(siteid == .x) %>%
        select(siteid, year, jaccard),
      y = "jaccard", x = "year",
      smoothing_method = "lm"))
plot_grid(plotlist = p_jaccard_low_jaccard, ncol = 3)

p_pop_low_jaccard <-
  map(site_desc_jaccard_slope[
    (length(site_desc_jaccard_slope)-20):length(site_desc_jaccard_slope)],
    ~plot_temporal_population(
      filtered_dataset$measurement %>%
        filter(siteid == .x)
      ) +
    theme(legend.position = "none")
    )
plot_grid(plotlist = p_pop_low_jaccard, ncol = 3)
#> Warning: Removed 36 rows containing missing values (position_stack).
#> Warning: Removed 132 rows containing missing values (position_stack).
#> Warning: Removed 225 rows containing missing values (position_stack).
#> Warning: Removed 20 rows containing missing values (position_stack).
#> Warning: Removed 25 rows containing missing values (position_stack).
#> Warning: Removed 18 rows containing missing values (position_stack).
#> Warning: Removed 9 rows containing missing values (position_stack).
#> Warning: Removed 12 rows containing missing values (position_stack).
#> Removed 12 rows containing missing values (position_stack).
#> Warning: Removed 10 rows containing missing values (position_stack).
#> Warning: Removed 102 rows containing missing values (position_stack).
#> Warning: Removed 12 rows containing missing values (position_stack).
#> Warning: Removed 45 rows containing missing values (position_stack).
#> Warning: Removed 100 rows containing missing values (position_stack).
#> Removed 100 rows containing missing values (position_stack).
#> Warning: Removed 25 rows containing missing values (position_stack).
#> Warning: Removed 24 rows containing missing values (position_stack).
#> Warning: Removed 60 rows containing missing values (position_stack).
#> Warning: Removed 30 rows containing missing values (position_stack).
#> Warning: Removed 105 rows containing missing values (position_stack).

p_turnover_low_jaccard_slope <- map(site_desc_jaccard_slope[(length(site_desc_jaccard_slope)-20):length(site_desc_jaccard_slope)],
  ~plot_community_data(
    analysis_dataset %>%
    filter(siteid == .x) %>%
    select(siteid, year, turnover),
    y = "turnover", x = "year",
  smoothing_method = "lm"))
plot_grid(plotlist = p_turnover_low_jaccard_slope, ncol = 3)

p_nestedness_low_jaccard_slope <- map(site_desc_jaccard_slope[(length(site_desc_jaccard_slope)-20):length(site_desc_jaccard_slope)],
  ~plot_community_data(
    analysis_dataset %>%
    filter(siteid == .x) %>%
    select(siteid, year, nestedness),
    y = "nestedness", x = "year",
  smoothing_method = "lm"))
plot_grid(plotlist = p_nestedness_low_jaccard_slope, ncol = 3)

p_appearance_low_jaccard_slope <- map(site_desc_jaccard_slope[(length(site_desc_jaccard_slope)-20):length(site_desc_jaccard_slope)],
  ~plot_community_data(
    analysis_dataset %>%
    filter(siteid == .x) %>%
    select(siteid, year, appearance),
    y = "appearance", x = "year",
  smoothing_method = "lm"))
plot_grid(plotlist = p_appearance_low_jaccard_slope, ncol = 3)

p_disappearance_low_jaccard_slope <- map(site_desc_jaccard_slope[(length(site_desc_jaccard_slope)-20):length(site_desc_jaccard_slope)],
  ~plot_community_data(
    analysis_dataset %>%
    filter(siteid == .x) %>%
    select(siteid, year, disappearance),
    y = "disappearance", x = "year",
  smoothing_method = "lm"))
plot_grid(plotlist = p_disappearance_low_jaccard_slope, ncol = 3)

4.1 Turnover and nestedness

stable_rich_slope <- slope_classification_station %>%
  filter(siteid %in% stable_rich_siteid) %>%
select(response, siteid, linear_slope) %>%
pivot_wider(names_from = "response", values_from = "linear_slope") %>%
mutate(
  turn_minus_nest = turnover - nestedness,
  turn_sup_0 = turnover > 0,
  nest_sup_0 = nestedness > 0,
  turn_sup_nest = turnover > nestedness,
  appearance_sup_disappearance = appearance > disappearance,
  jaccard_sup_0 = jaccard > 0,
  richness_sup_0 = species_nb > 0 
  )

4.1.1 Blowes et al. (2019)

ti <- table(stable_rich_slope$richness_sup_0, stable_rich_slope$turn_sup_0,
  dnn = list("richness > 0", "turnover > 0"))
ti
#>             turnover > 0
#> richness > 0 FALSE TRUE
#>        FALSE   915 1370
#>        TRUE   1100 1420
ti <- table(stable_rich_slope$richness_sup_0, stable_rich_slope$nest_sup_0,
  dnn = list("richness > 0", "nestedness > 0"))
ti
#>             nestedness > 0
#> richness > 0 FALSE TRUE
#>        FALSE   769 1516
#>        TRUE    670 1850
sum(stable_rich_slope$jaccard_sup_0)
#> [1] 505
ti <- table(stable_rich_slope$jaccard_sup_0, stable_rich_slope$turn_sup_nest,
  dnn = list("jaccard > 0", "turnover > nestedness"))
ti
#>            turnover > nestedness
#> jaccard > 0 FALSE TRUE
#>       FALSE  2059 2241
#>       TRUE    169  336
ti <- table(stable_rich_slope$jaccard_sup_0,
  stable_rich_slope$appearance_sup_disappearance,
  dnn = list("jaccard > 0", "appearance > disappearance"))
ti
#>            appearance > disappearance
#> jaccard > 0 FALSE TRUE
#>       FALSE  2032 2268
#>       TRUE    324  181
stable_rich_slope %>%
pivot_longer(-siteid, names_to = "response", values_to = "linear_slope") %>%
filter(response %in% c("jaccard_dis", "turnover", "nestedness", "turn_minus_nest")) %>%
select(response, siteid, linear_slope) %>%
ggplot(aes(x = response, y = linear_slope)) +
geom_boxplot()

5 Classification matrix

#debugonce(get_matrix_classification)
get_matrix_classification <- function(dataset = NULL, x = NULL, y = NULL) {
  ti <- dataset %>%
    filter(response %in% c(x, y)) %>%
    select(siteid, response, shape_class) %>%
    pivot_wider(names_from = "response", values_from = "shape_class")

  table(
    ti[[x]],
    ti[[y]]
  )
}
ti <- get_matrix_classification(
  dataset = slope_classification_station,
  x = "species_nb",
  y = "hillebrand"
  )

corrplot::corrplot(as.matrix(ti) / rowSums(ti))

5.1 Analysis

5.2 Reproducibility

Reproducibility receipt 
## datetime
Sys.time()
#> [1] "2022-04-18 18:59:32 CDT"

## repository
if(requireNamespace('git2r', quietly = TRUE)) {
  git2r::repository()
} else {
  c(
    system2("git", args = c("log", "--name-status", "-1"), stdout = TRUE),
    system2("git", args = c("remote", "-v"), stdout = TRUE)
  )
}
#> Local:    main L:/alain/RivFishTimeBiodiversityFacets
#> Remote:   main @ origin (https://github.com/alaindanet/RivFishTimeBiodiversityFacets.git)
#> Head:     [2326d9e] 2022-04-18: center hft_93 because with interactions, main coeffients are interpreted as when hft_93 = 0 which does not have any meaning ecologically in our dataset.

## session info
sessionInfo()
#> R version 4.1.2 (2021-11-01)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19042)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.1252 
#> [2] LC_CTYPE=English_United States.1252   
#> [3] LC_MONETARY=English_United States.1252
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.1252    
#> 
#> attached base packages:
#> [1] parallel  stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] ecocom_0.0.0.9000       factoextra_1.0.7        ade4_1.7-18            
#>  [4] broom.mixed_0.2.7       dotwhisker_0.7.4        skimr_2.1.3            
#>  [7] iNEXT_2.0.20            clValid_0.7             cluster_2.1.2          
#> [10] tclust_1.4-2            ggeffects_1.1.1         report_0.5.1           
#> [13] see_0.6.9               correlation_0.8.0       modelbased_0.7.2       
#> [16] effectsize_0.6.0.1      parameters_0.16.0       performance_0.8.0      
#> [19] bayestestR_0.11.5       datawizard_0.3.0        insight_0.16.0         
#> [22] easystats_0.4.3         glmmTMB_1.1.3           inlatools_0.0.1.9001   
#> [25] INLA_21.11.22           sp_1.4-6                foreach_1.5.2          
#> [28] Matrix_1.4-0            slider_0.2.2            vegan_2.5-7            
#> [31] lattice_0.20-45         permute_0.9-7           codyn_2.0.5            
#> [34] janitor_2.1.0           viridis_0.6.2           viridisLite_0.4.0      
#> [37] cowplot_1.1.1           terra_1.5-17            rnaturalearthdata_0.1.0
#> [40] rnaturalearth_0.1.0     sf_1.0-6                rmarkdown_2.11         
#> [43] scales_1.1.1            kableExtra_1.3.4        here_1.0.1             
#> [46] lubridate_1.8.0         magrittr_2.0.2          forcats_0.5.1          
#> [49] stringr_1.4.0           dplyr_1.0.7             purrr_0.3.4            
#> [52] readr_2.1.2             tidyr_1.2.0             tibble_3.1.6           
#> [55] ggplot2_3.3.5           tidyverse_1.3.1         tarchetypes_0.4.1      
#> [58] future.callr_0.7.0      future_1.23.0           targets_0.10.0         
#> [61] conflicted_1.1.0       
#> 
#> loaded via a namespace (and not attached):
#>   [1] utf8_1.2.2          ggstance_0.3.5      tidyselect_1.1.2   
#>   [4] lme4_1.1-28         grid_4.1.2          munsell_0.5.0      
#>   [7] base64url_1.4       codetools_0.2-18    units_0.8-0        
#>  [10] withr_2.5.0         colorspace_2.0-2    highr_0.9          
#>  [13] knitr_1.37          rstudioapi_0.13     ggsignif_0.6.3     
#>  [16] Rttf2pt1_1.3.10     listenv_0.8.0       labeling_0.4.2     
#>  [19] emmeans_1.7.2       git2r_0.29.0        repr_1.1.4         
#>  [22] farver_2.1.0        rprojroot_2.0.2     coda_0.19-4        
#>  [25] parallelly_1.30.0   vctrs_0.3.8         generics_0.1.2     
#>  [28] xfun_0.29           R6_2.5.1            cachem_1.0.6       
#>  [31] assertthat_0.2.1    nnet_7.3-17         gtable_0.3.0       
#>  [34] globals_0.14.0      processx_3.5.2      rlang_1.0.1        
#>  [37] systemfonts_1.0.3   splines_4.1.2       rstatix_0.7.0      
#>  [40] extrafontdb_1.0     TMB_1.7.22          broom_0.7.12       
#>  [43] abind_1.4-5         yaml_2.2.2          reshape2_1.4.4     
#>  [46] modelr_0.1.8        backports_1.4.1     extrafont_0.17     
#>  [49] tools_4.1.2         bookdown_0.24       ellipsis_0.3.2     
#>  [52] jquerylib_0.1.4     proxy_0.4-26        Rcpp_1.0.8         
#>  [55] plyr_1.8.6          base64enc_0.1-3     classInt_0.4-3     
#>  [58] ps_1.6.0            ggpubr_0.4.0        haven_2.4.3        
#>  [61] ggrepel_0.9.1       hrbrthemes_0.8.0    fs_1.5.2           
#>  [64] furrr_0.2.3         data.table_1.14.2   warp_0.2.0         
#>  [67] reprex_2.0.1        mvtnorm_1.1-3       hms_1.1.1          
#>  [70] evaluate_0.14       xtable_1.8-4        readxl_1.3.1       
#>  [73] gridExtra_2.3       compiler_4.1.2      KernSmooth_2.23-20 
#>  [76] crayon_1.5.0        minqa_1.2.4         htmltools_0.5.2    
#>  [79] mgcv_1.8-38         tzdb_0.2.0          DBI_1.1.2          
#>  [82] corrplot_0.92       dbplyr_2.1.1        MASS_7.3-55        
#>  [85] boot_1.3-28         car_3.0-12          cli_3.1.1          
#>  [88] igraph_1.2.11       pkgconfig_2.0.3     numDeriv_2016.8-1.1
#>  [91] xml2_1.3.3          svglite_2.1.0       bslib_0.3.1        
#>  [94] ggcorrplot_0.1.3    webshot_0.5.2       estimability_1.3   
#>  [97] rvest_1.0.2         snakecase_0.11.0    callr_3.7.0        
#> [100] digest_0.6.29       cellranger_1.1.0    gdtools_0.2.4      
#> [103] nloptr_2.0.0        lifecycle_1.0.1     nlme_3.1-155       
#> [106] jsonlite_1.8.0      carData_3.0-5       fansi_1.0.2        
#> [109] pillar_1.7.0        fastmap_1.1.0       httr_1.4.2         
#> [112] emo_0.0.0.9000      glue_1.6.1          iterators_1.0.14   
#> [115] class_7.3-20        stringi_1.7.6       sass_0.4.0         
#> [118] memoise_2.0.1       e1071_1.7-9